
###reduced version
use.formula = as.formula(paste('public.ingroup.preference ~ ',master.formula,sep = ''))
use.formula.reduced = as.formula(paste('public.ingroup.preference ~ ',master.formula.reduced,sep = ''))
use.vars = c('public.ingroup.preference',master.use.vars)

##UO
use.dat.uo = dat.uo[,use.vars]
use.dat.uo = use.dat.uo[complete.cases(use.dat.uo),]
##N per city for weighted regression below
city.N = tapply(use.dat.uo$income,use.dat.uo$index.city,length)
city.N = as.data.frame(city.N)
city.N$index.city = row.names(city.N)
use.dat.uo = merge(use.dat.uo,
                   city.N)

use.dat.uo$index.city.numeric = as.numeric(use.dat.uo$index.city)

reg.uo.pub.reduced = lm(use.formula.reduced, data = use.dat.uo)
clustered.se = cl(use.dat.uo, reg.uo.pub.reduced, use.dat.uo$index.city)
reg.uo.pub.reduced$se = clustered.se
reg.uo.pub.reduced.weighted = lm(use.formula.reduced, data = use.dat.uo,
                                  weights = city.N)
brl.se <- brl(use.formula.reduced, dataframe=use.dat.uo, 
              clusterName="index.city.numeric", method="brl")[,2]
reg.uo.pub.reduced.brl = reg.uo.pub.reduced
reg.uo.pub.reduced.brl$se = brl.se


reg.uo.pub = lm(use.formula, data = use.dat.uo)
clustered.se = cl(use.dat.uo, reg.uo.pub, use.dat.uo$index.city)
reg.uo.pub$se = clustered.se
reg.uo.pub.weighted = lm(use.formula, data = use.dat.uo,
                          weights = city.N)
brl.se <- brl(use.formula, dataframe=use.dat.uo, 
              clusterName="index.city.numeric", method="brl")[,2]
reg.uo.pub.brl = reg.uo.pub
reg.uo.pub.brl$se = brl.se

##Non UO
use.dat.sec = dat.non.uo[,use.vars]
use.dat.sec = use.dat.sec[complete.cases(use.dat.sec),]

##N per city for weighted regression below
city.N = tapply(use.dat.sec$income,use.dat.sec$index.city,length)
city.N = as.data.frame(city.N)
city.N$index.city = row.names(city.N)
use.dat.sec = merge(use.dat.sec,
                    city.N)

use.dat.sec$index.city.numeric = as.numeric(use.dat.sec$index.city)

reg.sec.pub.reduced = lm(use.formula.reduced, data = use.dat.sec)
clustered.se = cl(use.dat.sec, reg.sec.pub.reduced, use.dat.sec$index.city)
reg.sec.pub.reduced$se = clustered.se
reg.sec.pub.reduced.weighted = lm(use.formula.reduced, data = use.dat.sec,
                                   weights = city.N)
brl.se <- brl(use.formula.reduced, dataframe=use.dat.sec, 
              clusterName="index.city.numeric", method="brl")[,2]
reg.sec.pub.reduced.brl = reg.sec.pub.reduced
reg.sec.pub.reduced.brl$se = brl.se

reg.sec.pub = lm(use.formula, data = use.dat.sec)
clustered.se = cl(use.dat.sec, reg.sec.pub, use.dat.sec$index.city)
reg.sec.pub$se = clustered.se
reg.sec.pub.weighted = lm(use.formula, data = use.dat.sec,
                           weights = city.N)
brl.se <- brl(use.formula, dataframe=use.dat.sec, 
              clusterName="index.city.numeric", method="brl")[,2]
reg.sec.pub.brl = reg.sec.pub
reg.sec.pub.brl$se = brl.se


###############
outtable = apsrtable(reg.uo.pub.reduced,
                     reg.uo.pub,
                     reg.sec.pub.reduced,
                     reg.sec.pub, 
                     Sweave = T,
                     coef.names = coef.names,
                     model.names = c('UO','UO','Non UO','Non UO'),
                     notes = ''
                )
writeLines(
  outtable, 'output/appendix/Table_A6.tex')


###############

outtable = apsrtable(reg.uo.pub.reduced.weighted,
                     reg.uo.pub.weighted,
                     reg.sec.pub.reduced.weighted,
                     reg.sec.pub.weighted, 
                     Sweave = T,
                     coef.names = coef.names,
                     model.names = c('UO','UO','Non UO','Non UO'),
                     notes = ''
                      )
writeLines(
  outtable, 'output/appendix/Table_A15.tex')

outtable = apsrtable(reg.uo.pub.reduced.brl,
                     reg.uo.pub.brl,
                     reg.sec.pub.reduced.brl,
                     reg.sec.pub.brl, 
                     Sweave = T,
                     coef.names = coef.names,
                     model.names = c('UO','UO','Non UO','Non UO'),
                     notes = ''
)
writeLines(
  outtable, 'output/appendix/Table_A18.tex')



###zelig for effects of segregation

##UO
reg.uo = zelig(use.formula, data = use.dat.uo,
		model = 'ls',
		cite = F)

#expected values
outputs.uo = matrix(nrow = length(dissim.scale.uo), ncol = 3)
for(i in dissim.scale.uo){
	row.i = which(dissim.scale.uo == i)
	outx = setx(reg.uo, outgroup.dissim.yeshiva = i,  outgroup.yeshiva = uo.homo.value,
		jerusalem = 0, 
		demo1.sex = sex, 
    demo1.age = age, 
    demo1.ethnicity = ethnicity,
		demo2.left_right = ideology,
		income = income.value, 
    college = college.value, 
		nonjewish_pcnt = nonjewish_pcnt.value,
    immigrant = 0
    )
  sims = sim(reg.uo,outx, sims = n.sims)
  outs = sims$qi[[1]]
	outputs.uo[row.i,1] = mean(outs)
	outputs.uo[row.i,2] = quantile(outs,low.conf)
	outputs.uo[row.i,3] = quantile(outs,high.conf)
	}


###Non UO
reg.sec = zelig(use.formula, data = use.dat.sec,
		model = 'ls',
		cite = F)

		
outputs.sec = matrix(nrow = length(dissim.scale.sec), ncol = 3)
for(i in dissim.scale.sec){
	row.i = which(dissim.scale.sec == i)
	outx = setx(reg.sec, outgroup.dissim.yeshiva = i, outgroup.yeshiva = sec.homo.value,
		jerusalem = 0, 
		demo1.sex = sex, 
		demo1.age = age, 
		demo1.ethnicity = ethnicity,
		demo2.left_right = ideology,
		income = income.value, 
		college = college.value, 
		nonjewish_pcnt = nonjewish_pcnt.value,
		immigrant = 0
	)
	sims = sim(reg.sec,outx, sims = n.sims)
  out = summary(sims)$stats["Expected Values: E(Y|X)"]
	outs = sims$qi[[1]]
	outputs.sec[row.i,1] = mean(outs)
	outputs.sec[row.i,2] = quantile(outs,low.conf)
	outputs.sec[row.i,3] = quantile(outs,high.conf)
	}

	
###UO	
pdf('output/Figure_2_strategy_seg_uo.pdf',
    width = pdf.width,
    height = pdf.height)      
par(mar = c(4.5,3,.25,.25))
par(las = 1)  


plot(0,0,
	ylim = ylims,
	xlim = xlims.dissim.uo,
	type = 'n',
	xlab = 'Segregation',
  ylab = '',
  cex.axis = cex.axis,
  cex.lab = cex.lab)

abline(h = 0,
		col = 'gray',
		lty = 2)
lines(dissim.scale.uo,
	outputs.uo[,1],
	lwd = 4,
	col = 'black'
	)
lines(dissim.scale.uo,
	outputs.uo[,2],
	lty = 3,
	lwd = 3,
	col = 'black'
	)
lines(dissim.scale.uo,
	outputs.uo[,3],
	lty = 3,
	lwd = 3,
	col = 'black'
	)
	
dev.off()
	
##sec
pdf('output/Figure_2_strategy_seg_non_uo.pdf',
    width = pdf.width,
    height = pdf.height)      
#par(mar = c(4.5,5.5,.25,.25))
par(mar = c(3,3,.25,.25))

par(las = 1)  

plot(0,0,
	ylim = ylims,
	xlim = xlims.dissim.sec,
	type = 'n',
#	xlab = 'Segregation',
#	ylab = 'Ingroup Preference',
  xlab = '',
  ylab = '',
  cex.axis = cex.axis,
	cex.lab = cex.lab)
abline(h = 0,
		col = 'gray',
		lty = 2)
lines(dissim.scale.sec,
	outputs.sec[,1],
	lwd = 4,
	col = 'black'
	)
lines(dissim.scale.sec,
	outputs.sec[,2],
	lty = 3,
	lwd = 3,
	col = 'black'
	)
lines(dissim.scale.sec,
	outputs.sec[,3],
	lty = 3,
	lwd = 3,
	col = 'black'
	)
	
dev.off()


#####################
###zelig for effects of homogeneity


##UO
outputs.uo = matrix(nrow = length(homo.scale.uo), ncol = 3)
for(i in homo.scale.uo){
	row.i = which(homo.scale.uo == i)
	outx = setx(reg.uo,  outgroup.yeshiva = i, outgroup.dissim.yeshiva = mean.seg.value,
		jerusalem = 0, 
		demo1.sex = sex, 
		demo1.age = age, 
		demo1.ethnicity = ethnicity,
		demo2.left_right = ideology,
		income = income.value, 
		college = college.value, 
		nonjewish_pcnt = nonjewish_pcnt.value,
		immigrant = 0
	)
	sims = sim(reg.uo,outx, sims = n.sims)
	out = summary(sims)$stats["Expected Values: E(Y|X)"]
	outs = sims$qi[[1]]
	outputs.uo[row.i,1] = mean(outs)
	outputs.uo[row.i,2] = quantile(outs,low.conf)
	outputs.uo[row.i,3] = quantile(outs,high.conf)
	}


##Non UO	
outputs.sec = matrix(nrow = length(homo.scale.sec), ncol = 3)
for(i in homo.scale.sec){
	row.i = which(homo.scale.sec == i)
	outx = setx(reg.uo, outgroup.yeshiva = i, outgroup.dissim.yeshiva = mean.seg.value,
		jerusalem = 0, 
		demo1.sex = sex, 
		demo1.age = age, 
		demo1.ethnicity = ethnicity,
		demo2.left_right = ideology,
		income = income.value, 
		college = college.value, 
		nonjewish_pcnt = nonjewish_pcnt.value,
		immigrant = 0
	)
	sims = sim(reg.sec,outx, sims = n.sims)
	out = summary(sims)$stats["Expected Values: E(Y|X)"]
	outs = sims$qi[[1]]
	outputs.sec[row.i,1] = mean(outs)
	outputs.sec[row.i,2] = quantile(outs,low.conf)
	outputs.sec[row.i,3] = quantile(outs,high.conf)
	}

	
#UO	
pdf('output/Figure_2_strategy_proportion_uo.pdf',
    width = pdf.width,
    height = pdf.height)      
par(mar = c(4.5,3,.25,.25))

par(las = 1)  
plot(0,0,
	ylim = ylims,
	xlim = xlims.homo.uo,
	type = 'n',
	xlab = 'Outgroup Proportion',
  ylab = '',
  cex.axis = cex.axis,
	cex.lab = cex.lab)

abline(h = 0,
		col = 'gray',
		lty = 2)
lines(homo.scale.uo,
	outputs.uo[,1],
	lwd = 4,
	col = 'black'
	)
lines(homo.scale.uo,
	outputs.uo[,2],
	lty = 3,
	lwd = 3,
	col = 'black'
	)
lines(homo.scale.uo,
	outputs.uo[,3],
	lty = 3,
	lwd = 3,
	col = 'black'
	)

dev.off()
	
##sec
pdf('output/Figure_2_strategy_proportion_non_uo.pdf',
    width = pdf.width,
    height = pdf.height)      
par(mar = c(3,3,.25,.25))

par(las = 1)	
plot(0,0,
	ylim = ylims,
	xlim = xlims.homo.sec,
	type = 'n',
  ylab = ''	,
  xlab = '',
	cex.axis = cex.axis,
	cex.lab = cex.lab)
  
abline(h = 0,
		col = 'gray',
		lty = 2)


lines(homo.scale.sec,
	outputs.sec[,1],
	lwd = 4,
	col = 'black'
	)
lines(homo.scale.sec,
	outputs.sec[,2],
	lty = 3,
	lwd = 3,
	col = 'black'
	)
lines(homo.scale.sec,
	outputs.sec[,3],
	lty = 3,
	lwd = 3,
	col = 'black'
	)
	
dev.off()

